
data {
  int<lower=1> N; // number of observations
  int<lower=1> R; // number of suspect races
  int<lower=1> D; // number of counties
  int<lower=1> K; // number of departments
  
  int<lower=1,upper=R> r[N]; // race of suspect
  int<lower=1,upper=D> d[N]; // county where stop occurred
  int<lower=1,upper=K> k[N]; // department where stop occurred
  
  int<lower=1> n[N]; // number of stops
    int<lower=0> s[N]; // number of searches
    int<lower=0> h[N]; // number of successful searches (hits)
}

parameters {	
  // hyperparameters
  real<lower=0> sigma_t; // standard deviation for the normal the thresholds are drawn from. 
  
  // search thresholds
  vector[N] t_i;
  
  // parameters for signal distribution. 
  vector[R] phi_r;
  vector[D] phi_d;
  vector[K] phi_k;

  vector[R] delta_r; 
  vector[D] delta_d;
  vector[K] delta_k;
}

transformed parameters {
  vector[N] phi;
  vector[N] delta;
  vector<lower=0, upper=1>[N] search_rate;
  vector<lower=0, upper=1>[N] hit_rate;
  real successful_search_rate;
  real unsuccessful_search_rate;
  
  for (i in 1:N) {	
    // phi is the fraction of people of race r, d who are guilty (ie, carrying contraband)
    phi[i]    = inv_logit(phi_r[r[i]] + phi_d[d[i]] + phi_k[k[i]]);
    
    // delta is the center of the guilty distribution. 
    delta[i] = exp(delta_r[r[i]] + delta_d[d[i]] + delta_k[k[i]]);
    
    successful_search_rate = phi[i] * (1 - normal_cdf(t_i[i], delta[i], 1));
    unsuccessful_search_rate = (1 - phi[i]) * (1 - normal_cdf(t_i[i], 0, 1)); 
    search_rate[i] = (successful_search_rate + unsuccessful_search_rate);
    hit_rate[i] = successful_search_rate / search_rate[i];
   }
}

model {  
  // Draw threshold hyperparameters
  sigma_t ~ normal(0, 2);
  
  // Draw race parameters
  phi_r    ~ normal(-4, .5);
  delta_r ~ normal(0, .5);
  
  // Draw k parameters. 
  phi_k    ~ normal(0, .5);
  delta_k ~ normal(0, .5);
  
  // Draw department parameters
  phi_d    ~ normal(0, .1);    
  delta_d ~ normal(0, .1);    
  
  t_i ~ normal(0, sigma_t);
  
  s ~ binomial(n, search_rate);
  h ~ binomial(s, hit_rate);
  
}


generated quantities {}

